home *** CD-ROM | disk | FTP | other *** search
/ The 640 MEG Shareware Studio 2 / The 640 Meg Shareware Studio CD-ROM Volume II (Data Express)(1993).ISO / prog / yamp2.zip / VIRTOP.CPP < prev    next >
C/C++ Source or Header  |  1992-01-16  |  45KB  |  1,628 lines
  1.  
  2. /*
  3. YAMP - Yet Another Matrix Program
  4. Version: 1.2
  5. Author: Mark Von Tress, Ph.D.
  6. Date: 01/23/92
  7.  
  8. Copyright(c) Mark Von Tress 1992
  9. Portions of this code are (c) 1991 by Allen I. Holub and are used by
  10. permission
  11.  
  12. DISCLAIMER: THIS PROGRAM IS PROVIDED AS IS, WITHOUT ANY
  13. WARRANTY, EXPRESSED OR IMPLIED, INCLUDING BUT NOT LIMITED
  14. TO FITNESS FOR A PARTICULAR PURPOSE. THE AUTHOR DISCLAIMS
  15. ALL LIABILITY FOR DIRECT OR CONSEQUENTIAL DAMAGES RESULTING
  16. FROM USE OF THIS PROGRAM.
  17.  
  18. */
  19.  
  20. #include "virt.h"
  21. /////////////////// matrix functions ---- non-member functions
  22.  
  23. int Setwid(int w )
  24.   {
  25.      static int wid=6;
  26.      wid = ( w < 0 ) ? wid : w;
  27.      return wid;
  28.   }
  29. int Setdec ( int d )
  30.   {
  31.      static int dec = 3;
  32.      dec = ( d < 0 ) ? dec : d;
  33.      if ( d >= Getwid() ) dec = Getwid() - 1;
  34.      dec = (dec <=0 ) ? 0 : dec;
  35.      return dec;
  36.   }
  37. int Getwid( void )
  38.   {
  39.       static int wid=6;
  40.       wid = Setwid( -1 );
  41.       return wid;
  42.   }
  43. int Getdec(void )
  44.   {
  45.       static int dec=3;
  46.       dec = Setdec( -1 );
  47.       return dec;
  48.   }
  49.  
  50.  
  51. VMatrix &Reada( char *fid) /* this reads an ascii matrix */
  52.   {
  53.       unsigned int chh[MAXCHARS],*chc;
  54.       int i=0,j=0,x1=0,x2=0;
  55.       FILE *file;
  56.       char mname[MAXCHARS];
  57.       float x=0;
  58.  
  59.  
  60.       file = fopen( fid, "r");
  61.       if( !file ) Dispatch->Nrerror("READA: error while opening file");
  62.  
  63.       for( i = 0; i<80; i++) chh[i] = 0;
  64.  
  65.       i = 0;
  66.  
  67.       do { /* read until first carriage return */
  68.          j = fgetc( file );
  69.      chh[i] = j;
  70.          i++;
  71.       }    while( j != ((int) 0x0A) && !ferror(file) && (i<80)) ;
  72.  
  73.      if ( ferror(file) ) Dispatch->Nrerror("READA: error while reading matrix name");
  74.  
  75.      chh[i] = (int) '\0'; /* null terminate the string */
  76.  
  77.      if ( chh[0] != (int) 0x0A ) {
  78.       strcpy( mname,((const char *) chh ));
  79.       chc = chh;
  80.       for(j = 1; j<i-1; j++) strcat( mname, (( const char *) (++chc) ));
  81.      }
  82.      else
  83.      {
  84.       strcpy(mname," ");
  85.      }
  86.  
  87.       fscanf(file,"%d%d",&x1,&x2);
  88.  
  89.       VMatrix *mat =new VMatrix(mname,x1,x2);
  90.       for( i=1; i<=x1; i++){
  91.      for( j=1; j<=x2; j++){
  92.         fscanf(file,"%f",&x);
  93.         mat->M(i,j) = (double) x ;
  94.      }
  95.       }
  96.       fclose( file );
  97.       Dispatch->Push( *mat );
  98.       delete mat;
  99.       return Dispatch->ReturnMat();
  100.   }   /* reada */
  101.  
  102. VMatrix &Readb(char *fid) /* this reads an binary matrix */
  103.   {
  104.       int i=0,j=0,x1=0,x2=0;
  105.       FILE *file;
  106.       char name[MAXCHARS];
  107.       double x=0;
  108.  
  109.       file = fopen( fid, "rb");
  110.       if( !file ) Dispatch->Nrerror("READB: error while opening file");
  111.  
  112.       fread(&i, sizeof(int), 1, file);
  113.       fgets( name, (i+1), file);
  114.       fread(&x1, sizeof(int), 1, file);
  115.       fread(&x2, sizeof(int), 1, file);
  116.  
  117.       VMatrix *mat = new VMatrix( name,x1,x2);
  118.       for( i=1; i<=x1; i++){
  119.        for( j=1; j<=x2; j++){
  120.         fread(&x, sizeof(double), 1, file);
  121.         mat->M(i,j) = x ;
  122.      }
  123.       }
  124.       fclose( file );
  125.       Dispatch->Push( *mat );
  126.       delete mat;
  127.       return Dispatch->ReturnMat();
  128.   }   /* readb */
  129.  
  130. void Writea(char *fid, VMatrix &mat) /* write a matrix in an ascii file */
  131.   {
  132.      int i,j;
  133.      FILE *file;
  134.  
  135.      file = fopen( fid, "w" );
  136.      if ( !file ) Dispatch->Nrerror("WRITEA: unable to open file");
  137.  
  138.      mat.Garbage("Writea");
  139.  
  140.      fprintf(file,"%s\n", mat.StringAddress() );
  141.      fprintf(file,"%d %d \n", mat.r, mat.c );
  142.      for( i=1; i<= mat.r; i++){
  143.       for(j=1; j<= mat.c; j++)
  144.      fprintf(file,"%lf ",mat.m(i,j));
  145.       fprintf(file,"\n");
  146.      }
  147.      fclose(file);
  148.   } /* writea */
  149.  
  150. void heapsort( VMatrix &a )
  151. //  void svsort(n,v1,v2)
  152. //* sort uniform[0,1]'s in v1 and swap integers in
  153. //   v2 as uniforms are swapped */
  154. //    int n,v2[];
  155. //    double v1[];
  156.    {
  157.     int n=a.r;
  158.     int s0,s1,s2,s3,s4,s5,ts,s4m1,s5m1;
  159.     double t1;
  160.  
  161.     s0=s1= n; /* Shell-Metzger sort, PC-World, May 1986 */
  162.     s1 /= 2;
  163.     while(s1 > 0) {
  164.        s2=s0-s1;
  165.        for(s3=1; s3<=s2; s3++) {
  166.           s4=s3;
  167.           while (s4>0){
  168.             s5=s4+s1;
  169.             s4m1=s4;
  170.             s5m1=s5;
  171.             if ( a.m(s4m1,1) > a.m(s5m1,1) ){
  172.                t1=a.m(s4m1,1);
  173.                a.M(s4m1,1)=a.m(s5m1,1);
  174.                a.M(s5m1,1)=t1;
  175.                ts=a.m(s4m1,2);
  176.                a.M(s4m1,2)=a.m(s5m1,2);
  177.                a.M(s5m1,2)=ts;
  178.                s4=s4-s1;
  179.             }
  180.             else{
  181.                s4=0;
  182.             }
  183.           }
  184.        }
  185.        s1 /=2;
  186.     }
  187.    }
  188.  
  189.  
  190.  
  191. VMatrix &MSort(VMatrix &a, int col, int order)
  192.   {
  193.      a.Garbage();
  194.      VMatrix *sorted = new VMatrix("sorted(",a.r,a.c);
  195.      if( !order ){
  196.        // sort column col, then reorder other rows
  197.        VMatrix *temp = new VMatrix("temp",a.r,2);
  198.        for (int i=1; i<= a.r; i++ ) {
  199.      temp->M(i,1) = a.m(i,col);
  200.      temp->M(i,2) = (double)i;
  201.        }
  202.        heapsort( *temp );
  203.        for( i=1; i<=a.r; i++){
  204.       for( int j=1; j<=a.c; j++)
  205.            sorted->M( i,j ) = a.m( ((int)temp->m(i,2)), j );
  206.        }
  207.        delete temp;
  208.      }
  209.      else{
  210.        // sort row col, then reorder other columns
  211.        VMatrix *temp = new VMatrix("temp",a.c,2);
  212.        for (int i=1; i<= a.c; i++ ) {
  213.      temp->M(i,1) = a.m(col,i);
  214.      temp->M(i,2) = (double) i;
  215.        }
  216.        heapsort( *temp );
  217.        for( i=1; i<=a.c; i++){
  218.       for( int j=1; j<=a.r; j++)
  219.            sorted->M( j,i ) = a.m(j, ((int)temp->m(i,2)) );
  220.        }
  221.        delete temp;
  222.      }
  223.      strtype name = sorted->Getname(*sorted);
  224.      name = name+a.Getname(a)+")";
  225.      sorted->Nameit(name);
  226.      Dispatch->Push(*sorted);
  227.      delete sorted;
  228.      return Dispatch->ReturnMat();
  229.    }
  230.  
  231.  
  232.  
  233.  
  234.  
  235. VMatrix &Submat(VMatrix &a, int r, int c, int lr, int lc )
  236.    /* extract a submatrix of a into b*/
  237.   {
  238.       a.Garbage("Submat");
  239.  
  240.       if( (r > a.r ) || ( c > a.c) )
  241.     a.Nrerror( "SUBMAT: index out of range");
  242.  
  243.       VMatrix *b = new VMatrix( "b", r-lr+1, c-lc+1);
  244.       strtype temp = a.Getname(a);
  245.       b->Nameit( temp );
  246.       int r2 = r-lr+1;
  247.       int c2 = c-lc+1;
  248.       for(int i=1; i<=r2; i++){
  249.      for(int j=1; j<=c2; j++) b->M(i,j) = a.m(lr-1+i,lc-1+j);
  250.       }
  251.       Dispatch->Push(*b);
  252.       delete b;
  253.       return Dispatch->ReturnMat();
  254.   } /* submat */
  255.  
  256.  
  257. VMatrix &Ch(VMatrix &a,VMatrix &b ) /* horizontal concatination */
  258.   {
  259.      a.Garbage("Ch");
  260.      b.Garbage("Ch");
  261.      int adim = a.c;
  262.      int bdim = b.c;
  263.      int k = adim+bdim;
  264.  
  265.      if ( a.r != b.r)
  266.        a.Nrerror("CH: matrices have different number of rows");
  267.  
  268.      VMatrix *c = new VMatrix("(",a.r, k);
  269.      if( !c ) a.Nrerror("CH: failed to create temp storage");
  270.  
  271.      strtype mname = c->Getname(*c)+a.Getname(a)+"||"+b.Getname(b)+")";
  272.      c->Nameit( mname );
  273.  
  274.      for (   int i=1; i <= a.r; i++){
  275.      for(int j=1; j <= a.c; j++)
  276.         c->M(i,j) = a.m(i,j);
  277.      }
  278.      int ii;
  279.      for (        i=1, ii=1; i<=b.r; i++, ii++){
  280.      for( int j=1, jj=1; j<=b.c; j++, jj++) c->M(ii,jj+a.c) = b.m(i,j);
  281.      }
  282.      Dispatch->Push(*c);
  283.      delete c;
  284.      return Dispatch->ReturnMat();
  285.   } /* ch */
  286.  
  287. VMatrix &Cv(VMatrix &a,VMatrix &b ) /* vertical concatination */
  288.   {
  289.      a.Garbage("Cv");
  290.      b.Garbage("Cv");
  291.      int adim = a.r;
  292.      int bdim = b.r;
  293.      int k = adim + bdim;
  294.  
  295.      if ( a.c != b.c)
  296.        a.Nrerror("CV: matrices have different number of columns");
  297.  
  298.      VMatrix *c = new VMatrix("(",k, a.c);
  299.  
  300.      strtype mname = c->Getname(*c)+a.Getname(a)+"//"+b.Getname(b)+")";
  301.      c->Nameit( mname );
  302.  
  303.      for (    int i= 1; i <= a.r; i++){
  304.      for( int j= 1; j <= a.c; j++)
  305.         c->M(i,j) =  a.m(i,j);
  306.      }
  307.      int ii;
  308.      for (        i=1, ii=1;  i<=b.r; i++, ii++){
  309.      for( int j=1, jj=1;  j<=b.c; j++, jj++) c->M(ii+a.r,jj) = b.m(i,j);
  310.      }
  311.      Dispatch->Push(*c);
  312.      delete c;
  313.      return Dispatch->ReturnMat();
  314.   } /* cv */
  315.  
  316.  
  317. //////////////////  math funtions
  318.  
  319. VMatrix& Tran(VMatrix &ROp )
  320.   {
  321.     ROp.Garbage("Tran");
  322.     VMatrix *temp = new VMatrix("t",ROp.c,ROp.r);
  323.     strtype newname = temp->Getname(ROp);
  324.     for( int i=1; i<=ROp.r; i++){
  325.       for( int j=1; j<=ROp.c; j++) temp->M(j,i) = ROp.m(i,j);
  326.     }
  327.     newname = newname + "'";
  328.     temp->Nameit( newname ) ;
  329.     Dispatch->Push(*temp);
  330.     delete temp;
  331.     return Dispatch->ReturnMat();
  332.   }
  333. VMatrix& Mexp(VMatrix &ROp )
  334.   {
  335.     // take exponent of all elements
  336.     ROp.Garbage("Mexp");
  337.     VMatrix *temp = new VMatrix("exp(",ROp.r,ROp.c);
  338.     for( int i=1; i<=ROp.r; i++){
  339.       for( int j=1; j<=ROp.c; j++) temp->M(i,j) = exp(ROp.m(i,j));
  340.     }
  341.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  342.     temp->Nameit( newname ) ;
  343.     Dispatch->Push(*temp);
  344.     delete temp;
  345.     return Dispatch->ReturnMat();
  346.   }
  347. VMatrix& Mabs(VMatrix &ROp )
  348.   {
  349.     // take absolute values of all elements
  350.     ROp.Garbage("Mabs");
  351.     VMatrix *temp = new VMatrix("abs(",ROp.r,ROp.c);
  352.     for( int i=1; i<=ROp.r; i++){
  353.       for( int j=1; j<=ROp.c; j++) temp->M(i,j) = fabs(ROp.m(i,j));
  354.     }
  355.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  356.     temp->Nameit( newname ) ;
  357.     Dispatch->Push(*temp);
  358.     delete temp;
  359.     return Dispatch->ReturnMat();
  360.   }
  361. VMatrix& Mlog(VMatrix &ROp )
  362.   {
  363.     // take logarithm of all elements
  364.     ROp.Garbage("Mlog");
  365.     VMatrix *temp = new VMatrix("log(",ROp.r,ROp.c);
  366.     for( int i=1; i<=ROp.r; i++){
  367.       for( int j=1; j<=ROp.c; j++) {
  368.      double R = ROp.m(i,j);
  369.      temp->M(i,j) = ( R > 0.0 ) ? log( R ) :
  370.          Dispatch->Nrerror("Mlog: zero or negative argument to log");
  371.       }
  372.     }
  373.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  374.     temp->Nameit( newname ) ;
  375.     Dispatch->Push(*temp);
  376.     delete temp;
  377.     return Dispatch->ReturnMat();
  378.   }
  379. VMatrix& Msqrt(VMatrix &ROp )
  380. {
  381.     // take sqrt of all elements
  382.     ROp.Garbage("Msqrt");
  383.     VMatrix *temp = new VMatrix("sqrt(",ROp.r,ROp.c);
  384.     for( int i=1; i<=ROp.r; i++){
  385.       for( int j=1; j<=ROp.c; j++) {
  386.      double R = ROp.m(i,j);
  387.      temp->M(i,j) = ( R >= 0.0 ) ? sqrt( R ) :
  388.          Dispatch->Nrerror("Msqrt: zero or negative argument to sqrt");
  389.       }
  390.     }
  391.     strtype newname = temp->Getname(*temp) + ROp.Getname(ROp) + ")";
  392.     temp->Nameit( newname ) ;
  393.     Dispatch->Push(*temp);
  394.     delete temp;
  395.     return Dispatch->ReturnMat();
  396.   }
  397. double Trace( VMatrix &ROp )
  398.   {
  399.      ROp.Garbage("Trace");
  400.      double trace = 0;
  401.      if( ROp.r != ROp.c )
  402.        ROp.Nrerror("Trace: matrix not square in trace");
  403.      for( int i=1; i<=ROp.r; i++) trace += ROp.m(i,i);
  404.      return trace;
  405.   }
  406. VMatrix& Ident(int n )
  407.   {
  408.     VMatrix *temp = new VMatrix("I(",n,n);
  409.     for( int i=1; i<=n; i++) {
  410.        for( int j=1; j<=n; j++) temp->M(i,j) =(i == j) ? 1.0 : 0.0;
  411.     }
  412.  
  413.     char buffer[MAXCHARS]={'\0'};
  414.     gcvt( ((double) n), NDECS, buffer);
  415.     strtype string = buffer;
  416.     string = temp->Getname(*temp)+string+")";
  417.  
  418.     temp->Nameit( string );
  419.     Dispatch->Push(*temp);
  420.     delete temp;
  421.     return Dispatch->ReturnMat();
  422.   }
  423. VMatrix& Helm( int n )
  424.   {
  425.     VMatrix *temp = new VMatrix("Helm(",n,n);
  426.     for( int i=1; i<=n; i++) {
  427.        for( int j=1; j<=n; j++){
  428.      double d = 1.0/sqrt( (double) n );
  429.      double den = ( j>1 ) ? 1.0 / sqrt( (double) j*(j-1) ) : d;
  430.      if( j>1 && j<i ) d = 0;
  431.      if( j>1 && j==i) d =  - den * ( (double ) (j-1) );
  432.      if( j>1 && j>i ) d = den;
  433.       temp->M(i,j) = d;
  434.      }
  435.     }
  436.     char buffer[MAXCHARS]={'\0'};
  437.     gcvt( ((double) n), NDECS, buffer);
  438.     strtype string = buffer;
  439.     string = temp->Getname(*temp)+string+")";
  440.     temp->Nameit( string );
  441.     Dispatch->Push(*temp);
  442.     delete temp;
  443.     return Dispatch->ReturnMat();
  444.   }
  445. VMatrix& Fill(int r, int c, double d )
  446.   {
  447.     VMatrix *temp = new VMatrix("j",r,c);
  448.     for( int i=1; i<=r; i++) {
  449.     for( int j=1; j<= c; j++) temp->M(i,j) = d;
  450.     }
  451.     Dispatch->Push(*temp);
  452.     delete temp;
  453.     return Dispatch->ReturnMat();
  454.   }
  455. VMatrix& Index( int lower, int upper, int step)
  456.   {
  457.      if( step == 0 ) Dispatch->Nrerror("Index: step is zero");
  458.      int cnter = 0;
  459.      if( lower < upper ){
  460.        for( int i=lower; i<= upper; cnter++, i += step );
  461.      }
  462.      else{
  463.        for( int i=upper; i>= lower; cnter++, i += step );
  464.      }
  465.      if (cnter == 0)
  466.        Dispatch->Nrerror("Index: trying to allocate a matrix with zero rows");
  467.      VMatrix *temp = new VMatrix("Index",cnter, 1);
  468.  
  469.      cnter = 1;
  470.      if( lower < upper ){
  471.        for( int i=lower; i<= upper; i += step )
  472.           temp->M( cnter++, 1) = (double) i;
  473.      }
  474.      else{
  475.        for( int i=upper; i>= lower; i += step )
  476.           temp->M( cnter++, 1) = (double) i;
  477.      }
  478.  
  479.      Dispatch->Push( *temp );
  480.      delete temp;
  481.      return Dispatch->ReturnMat();
  482.  }
  483.  
  484. VMatrix& Sweep( int k1, int k2, VMatrix& ROp) /* sweep out a row */
  485.   {
  486.      long double eps=ZERO, d;
  487.      int i,j,k,n, it;
  488.  
  489.      ROp.Garbage("Sweep");
  490.      if (ROp.r != ROp.c ) ROp.Nrerror("SWEEP: matrix a is not square");
  491.  
  492.      n = ROp.r;
  493.      if (k2 < k1 ){ k=k1; k1=k2; k2=k;}
  494.  
  495.      for ( k=k1; k<=k2; k++){
  496.      if( fabs(ROp.m(k,k)) < eps )
  497.         for( it = 1; it<= n; it++){
  498.         ROp.M(it,k) = 0;
  499.         ROp.M(k,it) = 0;
  500.      }
  501.      else {
  502.          d = 1.0/ROp.m(k,k);
  503.          ROp.M(k,k) = d;
  504.          for (i=1; i<=n; i++) {if ( i != k )
  505.          ROp.M(i,k) = ROp.m(i,k)*(double) -d;}
  506.          for (j=1; j<=n; j++) {if ( j != k )
  507.          ROp.M(k,j) = ROp.m(k,j)*(double)  d;}
  508.          for (i=1; i<=n; i++) {
  509.          if( i != k ) {
  510.              for(j =1; j<=n; j++) {
  511.             if (j != k )
  512.               ROp.M(i,j) = ROp.m(i,j)+(ROp.m(i,k))*(ROp.m(k,j))/d;
  513.              }
  514.          }
  515.          }
  516.      }
  517.      }
  518.      strtype newname = "sweep(";
  519.      newname = newname + ROp.Getname(ROp) + ")";
  520.      ROp.Nameit( newname );
  521.      Dispatch->Push( ROp );
  522.      return  Dispatch->ReturnMat();
  523.     }/*sweep*/
  524.  
  525. VMatrix& Inv( VMatrix &ROp ) /*matrix inversion using sweep */
  526.   {
  527.      struct mat *b;
  528.      ROp.Garbage("Inv");
  529.      Dispatch->Inclevel();
  530.  
  531.      if( ROp.c != ROp.r ) ROp.Nrerror("INVERSE: matrix a not square");
  532.      strtype newname = "invs(";
  533.      newname = newname + ROp.Getname(ROp) + ")";
  534.      VMatrix *temp = new VMatrix("t",ROp.c,ROp.r);
  535.  
  536.      *temp = ROp;
  537.      *temp = Sweep(1, temp->r, *temp);
  538.  
  539.      temp->Nameit( newname );
  540.      Dispatch->Push( *temp );
  541.      delete temp;
  542.      return Dispatch->DecReturn();
  543.   } /* inverse */
  544.  
  545. VMatrix &Kron(VMatrix &a, VMatrix &b)
  546.   {
  547.      a.Garbage("Kron");
  548.      b.Garbage("Kron");
  549.  
  550.      int cr = a.r*b.r;
  551.      int cc = a.c*b.c;
  552.      VMatrix *c = new VMatrix( "(", cr, cc);
  553.      strtype mname = c->Getname(*c)+a.Getname(a)+"@"+b.Getname(b)+")";
  554.      c->Nameit( mname );
  555.  
  556.      for(int i=1; i<=a.r; i++){
  557.        for(int j=1; j<=a.c; j++){
  558.          for(int k=1; k<=b.r; k++ ){
  559.            for(int l=1; l<=b.c; l++){
  560.               c->M((i-1)*b.r+k,(j-1)*b.c+l) = a.m(i,j)*b.m(k,l);
  561.            }
  562.          }
  563.        }
  564.      }
  565.      Dispatch->Push(*c);
  566.      delete c;
  567.      return Dispatch->ReturnMat();
  568.   } /* kron */
  569.  
  570.  
  571. void Pivot(VMatrix &Data, int RefRow,
  572.        double &Determ,enum boolean &DetEqualsZero)
  573. {
  574.   DetEqualsZero = TTRUE;
  575.   int NewRow = RefRow;
  576.   while (DetEqualsZero && (NewRow < Data.r) ) {
  577.     NewRow++;
  578.     if (fabs(Data.m(NewRow, RefRow)) > ZERO) {
  579.       for(int i=1; i<=Data.r; i++){
  580.         double temp = Data.m(NewRow,i);
  581.         Data.M(NewRow,i)=Data.m(RefRow,i);
  582.         Data.M(RefRow,i) = temp;
  583.       }
  584.       DetEqualsZero = FFALSE;
  585.       Determ = -Determ;  /* row swap adjustment to det */
  586.     }
  587.   }
  588. }
  589. double Det( VMatrix &Datain)
  590. {
  591.   Datain.Garbage("Det");
  592.  
  593.   if( Datain.r == Datain.c ){
  594.     Dispatch->Inclevel();
  595.     int Dimen = Datain.r;
  596.     VMatrix *Data = new VMatrix("data",Dimen,Dimen);
  597.  
  598.     *Data= Datain;
  599.     double Determ = 1;
  600.     long double Coef;
  601.     int  Row, RefRow = 0;
  602.     enum boolean  DetEqualsZero = FFALSE;
  603.  
  604.     while (!(DetEqualsZero) && (RefRow < Dimen - 1) ) {
  605.       RefRow++;
  606.       if (fabs(Data->m(RefRow, RefRow)) < ZERO)
  607.          Pivot(*Data, RefRow, Determ, DetEqualsZero);
  608.       if (!(DetEqualsZero))
  609.       for (Row = RefRow + 1; Row <= Dimen; Row++)
  610.        if (fabs(Data->m(Row, RefRow)) > ZERO) {
  611.            Coef = -((long double)Data->m(Row, RefRow)) /
  612.                    ((long double)Data->m(RefRow, RefRow));
  613.            for( int i=1; i<=Dimen; i++)
  614.               Data->M(Row,i) = Data->m(Row,i) +
  615.                 (double)(Coef*((long double)Data->m(RefRow,i)));
  616.        }
  617.      Determ *= Data->m(RefRow, RefRow);
  618.     }
  619.  
  620.     Determ =( DetEqualsZero ) ? 0 : Determ * Data->m(Dimen, Dimen);
  621.     delete Data;
  622.     Dispatch->Declevel();
  623.     return Determ;
  624.   }
  625.   else Datain.Nrerror("Det: Matrix is not square\n");
  626. }
  627.  
  628. //--------------- cholesky decomposition ----------------------
  629. // ROp is symmetric, returns upper triangular S where ROp = S'S
  630. //-------------------------------------------------------------
  631. #define EPS 1.0e-7
  632. VMatrix& Cholesky(VMatrix& ROp)
  633. {
  634.    int nr = ROp.r;
  635.  
  636.    ROp.Garbage("Cholesky");
  637.    for( int i=1; i<=nr; i++){
  638.      for( int j=1; j< i; j++)
  639.        if( fabs( ROp.m(i,j) - ROp.m(j,i) ) > EPS )
  640.          Dispatch->Nrerror("Cholesky: Input matrix is not symmetric");
  641.    }
  642.    VMatrix  *temp = new VMatrix();
  643.    *temp = Fill(nr,nr,0.0);
  644.    for ( i=1; i<=nr; i++){
  645.       for (int j = i; j<=nr; j++){
  646.          double sum = 0.0;
  647.          for (int k=1 ; k < i; k++) {
  648.             sum += temp->m(k,i) * temp->m(k,j);
  649.          }
  650.          if (i==j) temp->M(i,j) = ( ROp.m(i,j) < sum ) ?
  651.                    Dispatch->Nrerror("Cholesky: negative pivot") :
  652.                    sqrt(ROp.m(i,j) - sum);
  653.          else temp->M(i,j) = ( fabs(temp->m(i,i)) < EPS ) ?
  654.                   Dispatch->Nrerror("Cholesky: zero or negative pivot") :
  655.                   (ROp.m(i,j) - sum) / temp->m(i,i);
  656.       }
  657.    }
  658.    strtype newname =  "half(";
  659.    newname = newname + ROp.Getname(ROp) + ")";
  660.    temp->Nameit( newname );
  661.    Dispatch->Push(*temp);
  662.    delete temp;
  663.    return Dispatch->ReturnMat();
  664. }
  665.  
  666.  
  667. void QR( VMatrix& ROp, VMatrix& Q, VMatrix& R, boolean makeq)
  668.   {
  669.     ROp.Garbage("QR");
  670.     Q.Garbage("QR");
  671.     R.Garbage("QR");
  672.  
  673.     Dispatch->Inclevel();
  674.     int r=ROp.r;
  675.     int c=ROp.c;
  676.  
  677.     Q = ROp;
  678.     R = Fill(c,c,0.0);
  679.     double s;
  680.  
  681.     for(int j=1; j<=c; j++){
  682.        double sigma = 0.0;
  683.        for( int i=j; i<=r; i++) sigma += Q.m(i,j)*Q.m(i,j);
  684.        if ( fabs( sigma) <= 1.0e-10 ) Dispatch->Nrerror("QR: singular X");
  685.        R.M(j,j) = s = ( Q.m(j,j) < 0 ) ? sqrt( sigma ) : -sqrt( sigma );
  686.        double beta = 1.0/(s*Q.m(j,j)-sigma);
  687.        Q.M(j,j) = Q.m(j,j) - s;
  688.        for(int k = j+1; k<=c; k++){
  689.           double sum = 0.0;
  690.           for( int i=j; i<=r; i++) sum += Q.m(i,j)*Q.m(i,k);
  691.           sum *= beta;
  692.           for( i=j; i<=r; i++) Q.M(i,k) = Q.m(i,k) + Q.m(i,j)*sum;
  693.           R.M(j,k) = Q.m(j,k);
  694.        }
  695.   }
  696.  
  697.   if ( makeq ) Q = ROp*Ginv(R);
  698.  
  699.   Q.Nameit("Q");
  700.   R.Nameit("R");
  701.   Dispatch->Declevel();
  702. }
  703.  
  704.  
  705.  
  706.  
  707. #undef EPS
  708. //------------------- Singular Valued Decomposition ----------
  709. // from EISPACK SVD
  710. //------------------------------------------------------------
  711.  
  712. void  svdtemp(VMatrix &a, VMatrix &w, boolean matu, VMatrix &u,
  713.           boolean matv, VMatrix &v, int &ierr, VMatrix &rv1)
  714. {
  715.     Dispatch->Inclevel();
  716.     a.Garbage("SVD");
  717.     w.Garbage("SVD");
  718.     u.Garbage("SVD");
  719.     v.Garbage("SVD");
  720.     rv1.Garbage("SVD");
  721.  
  722.         int m  = a.r;
  723.         int n  = a.c;
  724.  
  725.        int i,j,k,l,ii,i1,kk,k1,ll,l1,its,mn;
  726. //      VMatrix a("a",m,n),w("w",n,1),
  727. //              u("u",m,m),v("v",n,n),rv1("rv1",n,1);
  728.        double c,f,g,h,s,x,y,z,eps,scale,machep,fss;
  729. //      boolean matu,matv;
  730.  
  731. //       ********** Machep is a machine dependent parameter specifying
  732. //              the relative precision of floatin point aritmetic.
  733. //
  734. //       **************
  735.  
  736.        machep = pow(16.0,(-5.0));
  737.  
  738.        ierr = 0;
  739.  
  740.        u = a;
  741. //     ********householder reduction to bi diagonal form ********
  742.        g=0.0;
  743.        scale=0.0;
  744.        x=0.0;
  745.  
  746.        for( i=1; i<=n; i++){
  747.               l=i+1;
  748.               rv1.M(i,1)=scale*g;
  749.               g=0.0;
  750.               s=0.0;
  751.               scale=0.0;
  752.               if (i > m) goto l210;
  753.  
  754.               for (k=i; k<=m; k++) scale += fabs(u.m(k,i));
  755.  
  756.               if (scale == 0.0) goto l210;
  757.  
  758.               for( k=i; k <= m; k++){
  759.                      u.M(k,i)=u.m(k,i)/scale;
  760.                      s += u.m(k,i)*u.m(k,i);
  761.               }
  762.  
  763.               f=u.m(i,i);
  764.               g=-(( f >= 0.0 ) ? fabs( sqrt(s) ) : -fabs(sqrt(s)));
  765.               h=f*g-s ;
  766.               u.M(i,i)=f-g;
  767.               if (i==n) goto l190;
  768.  
  769.               for( j=l; j<=n; j++){
  770.                      s=0.0;
  771.                      for(k=i; k<=m; k++) s += u.m(k,i)*u.m(k,j);
  772.                      f=s/h;
  773.                      for(k=i; k<=m; k++) u.M(k,j)=u.m(k,j)+f*u.m(k,i);
  774.               }
  775.  l190:              for( k=i; k<=m; k++) u.M(k,i)=scale*u.m(k,i);
  776.  
  777.  l210:              w.M(i,1)=scale*g;
  778.               g=0.0;
  779.               s=0.0;
  780.               scale=0.0;
  781.               if(i>m || i==n) goto l290;
  782.  
  783.               for( k=l; k<=n; k++) scale += fabs(u.m(i,k));
  784.  
  785.               if (scale == 0.0) goto l290;
  786.  
  787.               for(k=l; k<=n; k++){
  788.                      u.M(i,k) =u.m(i,k)/scale;
  789.                      s += u.m(i,k)*u.m(i,k);
  790.               }
  791.               f=u.m(i,l);
  792.               g=-( (f >= 0.0) ? fabs(sqrt(s)): -fabs(sqrt(s)) );
  793.               h=f*g-s;
  794.               u.M(i,l)=f-g;
  795.  
  796.               for( k=l; k<=n; k++) rv1.M(k,1) = u.m(i,k)/h;
  797.  
  798.               if(i == m) goto l270;
  799.  
  800.               for( j=l; j<=m; j++){
  801.                      s=0.0;
  802.                      for(k=l; k<=n; k++) s+= u.m(j,k)*u.m(i,k);
  803.                      for(k=l; k<=n; k++) u.M(j,k)=u.m(j,k)+s*rv1.m(k,1);
  804.               }
  805.  
  806.  l270:              for( k=l; k<=n; k++ ) u.M(i,k)=scale*u.m(i,k);
  807.  
  808.  l290:              x = ( x > fabs(w.m(i,1))+fabs(rv1.m(i,1))) ? x :
  809.                                    fabs(w.m(i,1))+fabs(rv1.m(i,1)) ;
  810.        }
  811.  
  812. //       ********** accumulation of right-hand transformations ********
  813.  
  814.        if (!matv) goto l410;
  815.  
  816.        for( ii=1; ii<=n; ii++){
  817.               i=n+1-ii;
  818.               if(i == n) goto l390;
  819.               if(g == 0.0) goto l360;
  820.  
  821. //               double division avoids possible underflow
  822.               for(j=l; j<=n; j++) v.M(j,i) =(u.m(i,j)/u.m(i,l))/g;
  823.  
  824.               for( j=l; j<=n; j++){
  825.                      s=0.0;
  826.                      for(k=l; k<=n; k++) s+=u.m(i,k)*v.m(k,j);
  827.                      for(k=l; k<=n; k++) v.M(k,j)=v.m(k,j)+s*v.m(k,i);
  828.               }
  829.  
  830.  l360:              for( j=l; j<=n; j++){
  831.                      v.M(i,j)=0.0;
  832.                      v.M(j,i)=0.0;
  833.               }
  834.  
  835.  l390:              v.M(i,i)=1.0;
  836.               g=rv1.m(i,1);
  837.               l=i;
  838.         }
  839. //       *************accumulation of left-hand transformations
  840.  l410:       if (!matu) goto l510;
  841. //       ************* for i= min(m,n) step -1 until 1 do ******
  842.        mn=n;
  843.        if (m<n)mn=m;
  844.  
  845.        for(ii=1; ii<=mn; ii++){
  846.               i=mn+1-ii;
  847.           l=i+1;
  848.               g=w.m(i,1);
  849.               if(i==n) goto l430;
  850.  
  851.               for(j=l; j<=n; j++) u.M(i,j)=0.0;
  852.  
  853.  l430:              if (g == 0.0) goto l475;
  854.               if (i==mn) goto l460;
  855.  
  856.               for(j=l; j<=n; j++){
  857.                      s=0.0;
  858.  
  859.              for(k=l;k<=m; k++) s+=u.m(k,i)*u.m(k,j);
  860. //                       double division avoids possible underflow
  861.                      f=(s/u.m(i,i))/g;
  862.                      for(k=i; k<=m; k++) u.M(k,j)=u.m(k,j)+f*u.m(k,i);
  863.                 }
  864.  
  865.  l460:              for(j=i; j<=m; j++) u.M(j,i)=u.m(j,i)/g;
  866.  
  867.               goto l490;
  868.  
  869.  l475:              for(j=i; j<=m; j++) u.M(j,i)=0.0;
  870.  
  871.  l490:              u.M(i,i)=u.m(i,i)+1.0;
  872.         }
  873.  
  874. //       ********** diagonalization of the bidiagonal form ***********
  875.  
  876.  l510:       eps=machep*x;
  877. //       ********** for k=n step -1 until 1 do ***********************
  878.         for(kk=1; kk<=n; kk++){
  879.               k1=n-kk;
  880.               k=k1+1;
  881.               its=0;
  882. //       ********** test for splitting .**********
  883. //                 for l=k step -1 until 1 do
  884.  l520:              for( ll=1; ll <= k; ll++){
  885.                      l1=k-ll;
  886.                      l=l1+1;
  887.                      if (fabs(rv1.m(l,1))<=eps) goto l565;
  888. //                       rv1(1) is always zero, so there is no exit
  889. //                     through the bottom of the loop *********
  890.                      if (fabs(w.m(l1,1))<=eps) goto l540;
  891.         }
  892. //       ************* cancellation of rv1(l) if l greater than 1******
  893.  l540:              c=0.0;
  894.               s=1.0;
  895.  
  896.               for( i=l; i<=k; i++){
  897.                      f=s*rv1.m(i,1);
  898.                      rv1.M(i,1)=c*rv1.m(i,1);
  899.                      if (fabs(f) <= eps) goto l565;
  900.                      g=w.m(i,1);
  901.                      h=sqrt(f*f+g*g);
  902.                      w.M(i,1)=h;
  903.                      c=g/h;
  904.                      s=-f/h;
  905.                      if(matu){ //go to 560
  906.                              for( j=1; j<=m; j++){
  907.                     y=u.m(j,l1);
  908.                                     z=u.m(j,i);
  909.                                     u.M(j,l1)=y*c +z*s;
  910.                                     u.M(j,i)=-y*s +z*c;
  911.                               }
  912.                         }
  913.                 }
  914. //       ************** test for convergence ********************
  915.  l565:              z=w.m(k,1);
  916.               if( l==k) goto l650;
  917. //       *************** if shift from bottom 2 by 2 minor *******
  918.               if (its == 30) goto l1000;
  919.           its=its+1;
  920.               x=w.m(l,1);
  921.               y=w.m(k1,1);
  922.               g=rv1.m(k1,1);
  923.               h=rv1.m(k,1);
  924.               f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
  925.               g=sqrt(f*f+1.0);
  926.               fss = ( f >= 0.0)? fabs(g):-fabs(g);
  927.               f=((x-z)*(x+z)+h*(y/(f+fss)-h))/x;
  928. //       **************** next qr transformation ***************
  929.               c=1.0;
  930.               s=1.0;
  931. //
  932.               for(i1=l; i1<=k1; i1++){
  933.                      i=i1+1;
  934.                      g=rv1.m(i,1);
  935.                      y=w.m(i,1);
  936.                      h=s*g;
  937.                      g=c*g;
  938.                      z=sqrt(f*f+h*h);
  939.              rv1.M(i1,1)=z;
  940.                      c=f/z;
  941.                      s=h/z;
  942.                      f=x*c+g*s;
  943.              g=-x*s+g*c;
  944.                      h=y*s;
  945.                      y=y*c;
  946.                      if (matv){ // goto l575;
  947.                              for( j=1; j<=n; j++){
  948.                                     x=v.m(j,i1);
  949.                                     z=v.m(j,i);
  950.                                     v.M(j,i1)=x*c+z*s;
  951.                                     v.M(j,i)=-x*s+z*c;
  952.                               }
  953.                         }
  954.                      z=sqrt(f*f+h*h);
  955.              w.M(i1,1)=z;
  956. //       ************ rotation can be arbitrary if z is zero *******
  957.                      if (z != 0.0){ // goto l580;
  958.                              c=f/z;
  959.                              s=h/z;
  960.                         }
  961. // l580:
  962.                      f=c*g +s*y;
  963.                      x=-s*g +c*y ;
  964.                      if(matu){ //goto l600;
  965.                               for( j=1; j<=m; j++){
  966.                                     y=u.m(j,i1);
  967.                     z=u.m(j,i);
  968.                                     u.M(j,i1)=y*c+z*s;
  969.                                     u.M(j,i)=-y*s+z*c;
  970.                                 }
  971.                         }
  972.  
  973.  l600:              x=x;} //continue
  974.  
  975.               rv1.M(l,1)=0.0;
  976.               rv1.M(k,1)=f;
  977.               w.M(k,1)=x;
  978.               goto l520;
  979. //       ************** convergence **************
  980.  l650:              if(z >= 0.0) goto l700;
  981. //       *************  w(k) is made non-negative *********
  982.               w.M(k,1)=-z;
  983.               if (!matv) goto l700;
  984.  
  985.               for( j=1; j<=n; j++) v.M(j,k)=-v.m(j,k);
  986.  
  987.  l700:  x=x;}
  988.  
  989.  
  990.        ierr = 0;
  991.        Dispatch->Declevel();
  992.        return;
  993. //        ***************** set error -- no convergence to a
  994. //              singular value after 30 iterations
  995.  l1000: ierr=k;
  996.  Dispatch->Declevel();
  997.  return;
  998. }
  999.  
  1000. int Svd( VMatrix &A, VMatrix &U, VMatrix &S, VMatrix &V,
  1001.           boolean makeu, boolean makev)
  1002.  {
  1003.  
  1004.     Dispatch->Inclevel();
  1005.  
  1006.     A.Garbage("SVD");
  1007.     VMatrix aa, uu, ss, vv, rv1;
  1008.     int ierr = 0;
  1009.  
  1010.     if( A.r < A.c ) aa = Tran(A);
  1011.     else aa = A;
  1012.  
  1013.     uu = Fill( aa.r, aa.r, 0.0);
  1014.     vv = Fill( aa.c, aa.c, 0.0);
  1015.     ss = Fill( aa.c, 1, 0.0);
  1016.     rv1= Fill( aa.c, 1, 0.0);
  1017.  
  1018.  
  1019.     svdtemp( aa, ss, makeu, uu, makev, vv, ierr, rv1);
  1020.  
  1021.     if( A.r < A.c ) {
  1022.        if( makeu ) U = vv;
  1023.        if( makev ) V = uu;
  1024.     } else {
  1025.        if( makeu ) U = uu;
  1026.        if( makev ) V = vv;
  1027.     }
  1028.     S = ss;
  1029.     Dispatch->Declevel();
  1030.     return ierr;
  1031.   }
  1032. //---------------- end svd ------------------------------------
  1033.  
  1034. VMatrix& Ginv( VMatrix& ROp )
  1035.   {
  1036.     ROp.Garbage("Ginv");
  1037.     Dispatch->Inclevel();
  1038.  
  1039.     VMatrix *u=new VMatrix;
  1040.     VMatrix *s=new VMatrix;
  1041.     VMatrix *v=new VMatrix;
  1042.     VMatrix *g=new VMatrix;
  1043.  
  1044.     Svd( ROp, *u, *s, *v);
  1045.     for( int i=1; i<=s->r; i++){
  1046.        double t=s->m(i,1);
  1047.        s->M(i,1) = ( fabs(t) <= 1.0e-10) ? 0.0 : 1.0/t;
  1048.     }
  1049.     *g = *v *Diag(*s)*Tran(*u);
  1050.     strtype newname = "ginv(";
  1051.     newname = newname + ROp.Getname(ROp) + ")";
  1052.     g->Nameit( newname );
  1053.     Dispatch->Push( *g );
  1054.     delete u;
  1055.     delete s;
  1056.     delete v;
  1057.     delete g;
  1058.     return Dispatch->DecReturn();
  1059.  }
  1060.  
  1061. //-------------------------- fft ------------------------------
  1062. // de Boor (1980) SIAM J SCI. STAT. COMPUT. V1 no 1, pp 173-178
  1063. // and NewMat by R. G. Davies a C++ matrix Package
  1064. //-------------------------------------------------------------
  1065.  
  1066. static void cossin(int n, int d, double& c, double& s)
  1067. // calculate cos(twopi*n/d) and sin(twopi*n/d)
  1068. // minimise roundoff error
  1069. {
  1070.    long n4 = n * 4;
  1071.    int sector = (int)floor( (double)n4 / (double)d + 0.5 );
  1072.    n4 -= sector * d;
  1073.    if (sector < 0) sector = 3 - (3 - sector) % 4; else sector %= 4;
  1074.    double ratio = 1.5707963267948966192 * (double)n4 / (double)d;
  1075.  
  1076.    switch (sector)
  1077.    {
  1078.    case 0: c =  cos(ratio);
  1079.            s =  sin(ratio);
  1080.            break;
  1081.    case 1: c = -sin(ratio);
  1082.            s =  cos(ratio);
  1083.            break;
  1084.    case 2: c = -cos(ratio);
  1085.            s = -sin(ratio);
  1086.            break;
  1087.    case 3: c =  sin(ratio);
  1088.            s = -cos(ratio);
  1089.            break;
  1090.    }
  1091. }
  1092.  
  1093. static void fftstep(VMatrix &AB, VMatrix &XY,
  1094.        int after, int now, int before)
  1095. {
  1096.  
  1097.    int gamma = after * before;
  1098.    int delta = now * after;
  1099.  
  1100.    double r_arg = 1.0, i_arg = 0.0;
  1101.  
  1102.    int x = 1;
  1103.    int m = AB.r - gamma;
  1104.  
  1105.    for (int j = 0; j < now; j++){
  1106.       int a  = 1;
  1107.       int x1 = x;
  1108.       x += after;
  1109.       for (int ia = 0; ia < after; ia++){
  1110.          // generate sins & cosines explicitly rather than iteratively
  1111.          // for more accuracy; but slower
  1112.          cossin(-(j*after+ia), delta, r_arg, i_arg);
  1113.  
  1114.          int a1 =  a++;
  1115.          int x2 = x1++;
  1116.          if (now==2){
  1117.             int ib = before;
  1118.             while (ib--){
  1119.                 int a2 = m + a1;
  1120.                 a1 += after;
  1121.                 double r_value = AB.m(a2,1);
  1122.                 double i_value = AB.m(a2,2);
  1123.                 XY.M(x2,1) = r_value*r_arg - i_value*i_arg + AB.m(a2-gamma,1);
  1124.                 XY.M(x2,2) = r_value*i_arg + i_value*r_arg + AB.m(a2-gamma,2);
  1125.                 x2 += delta;
  1126.             }
  1127.          }
  1128.          else {
  1129.             int ib = before;
  1130.             while (ib--){
  1131.                int a2 = m + a1;
  1132.                a1 += after;
  1133.                double r_value = AB.m(a2,1);
  1134.                double i_value = AB.m(a2,2);
  1135.                int in = now-1;
  1136.                while (in--){
  1137.                   a2 -= gamma;
  1138.                   double temp = r_value;
  1139.                   r_value = r_value*r_arg - i_value*i_arg + AB.m(a2,1);
  1140.                   i_value = temp*i_arg    + i_value*r_arg + AB.m(a2,2);
  1141.                }
  1142.                XY.M(x2,1) = r_value;
  1143.                XY.M(x2,2) = i_value;
  1144.                x2 += delta;
  1145.             }
  1146.          }
  1147.  
  1148.       }
  1149.    }
  1150. }
  1151.  
  1152.  
  1153. VMatrix& Fft( VMatrix &ROp, boolean INZEE )
  1154. // if INZEE == TTRUE  then calculate fft
  1155. // if INZEE == FFALSE then calculate inverse fft
  1156. {
  1157.    ROp.Garbage("FFT");
  1158.    if( ROp.c > 2 ) Dispatch->Nrerror("FFT: Input has more than two columns");
  1159.    int n = ROp.r;                     // length of arrays
  1160.    const int nextmx = 26;
  1161.    int prime[26] = { 2,3,5,7,11,13,17,19,23,29,31,37,
  1162.                     41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83,
  1163.                     89, 97, 101 };
  1164.    int after = 1;
  1165.    int before = n;
  1166.    int next = 0;
  1167.    boolean inzee = TTRUE;
  1168.  
  1169.    Dispatch->Inclevel();
  1170.  
  1171.    VMatrix *AB = new VMatrix("AB",n,2);
  1172.    VMatrix *XY = new VMatrix("XY",n,2);
  1173.  
  1174.    if( ROp.c == 1 ){
  1175.      for( int i=1; i<=n; i++){
  1176.         AB->M(i,1) = ROp.m(i,1);
  1177.         AB->M(i,2) = 0.0;
  1178.      }
  1179.    } else *AB = ROp;
  1180.    *XY = Fill(n,2,0.0);
  1181.  
  1182.    if( INZEE == FFALSE ){
  1183.      // take complex conjugate for ifft
  1184.      for( int i=1; i<=n; i++ ) AB->M(i,2) = -AB->m(i,2);
  1185.    }
  1186.  
  1187.    do
  1188.    {
  1189.       int now, b1;
  1190.       for (;;){
  1191.       if (next < nextmx) now = prime[next];
  1192.          b1 = before / now;
  1193.          if (b1 * now == before) break;
  1194.          next++;
  1195.          now += 2;
  1196.       }
  1197.       before = b1;
  1198.  
  1199.       if (inzee) fftstep(*AB, *XY, after, now, before);
  1200.          else    fftstep(*XY, *AB, after, now, before);
  1201.  
  1202.       inzee = (inzee == TTRUE) ? FFALSE : TTRUE;
  1203.       after *= now;
  1204.    }
  1205.    while (before != 1);
  1206.  
  1207.    if ( inzee == TTRUE ) *XY = *AB;
  1208.    // divide by n for ifft
  1209.    if ( INZEE == FFALSE) *XY = (*XY)/((double) n);
  1210.  
  1211.    strtype newname = ( INZEE == TTRUE ) ? "fft( " : "ifft(";
  1212.    newname = newname + ROp.Getname(ROp) + ")";
  1213.    XY->Nameit( newname );
  1214.  
  1215.    Dispatch->Push( *XY );
  1216.    delete XY;
  1217.    delete AB;
  1218.    return Dispatch->DecReturn();
  1219. }
  1220.  
  1221.  
  1222. ////////////////////  reshaping functions
  1223.  
  1224. VMatrix& Vec( VMatrix& ROp )
  1225.   {  // send columns of ROp to a vector
  1226.      ROp.Garbage("Vec");
  1227.  
  1228.      long lrc = ((long) ROp.r)*((long) ROp.c);
  1229.      if (lrc >= 32768 || lrc < 1)
  1230.         Dispatch->Nrerror("Vec: vec produces invalid indices");
  1231.      int r = ROp.r*ROp.c;
  1232.      int c = 1;
  1233.      VMatrix *temp = new VMatrix("t",r,c);
  1234.      int cnter=1;
  1235.      for( int j=1; j<=ROp.c; j++)
  1236.        for( int i=1; i<=ROp.r; i++) temp->M(cnter++,1) = ROp.m(i,j);
  1237.  
  1238.      strtype newname = "vec(";
  1239.      newname = newname + ROp.Getname(ROp) + ")";
  1240.      temp->Nameit( newname );
  1241.      Dispatch->Push( *temp );
  1242.      delete temp;
  1243.      return Dispatch->ReturnMat();
  1244.   }
  1245.  
  1246. VMatrix& Vecdiag( VMatrix& ROp )
  1247.   {  // Extract the diagonal from a square matrix and put in a vector
  1248.      ROp.Garbage("Vecdiag");
  1249.      if( ROp.r != ROp.c )
  1250.        Dispatch->Nrerror("Vecdiag: ROp is not square");
  1251.  
  1252.      VMatrix *temp = new VMatrix( "vdiag(", ROp.r, 1);
  1253.      for( int i=1; i<=ROp.r; i++) temp->M(i,1) = ROp.m(i,i);
  1254.  
  1255.      strtype newname = "vdiag(";
  1256.      newname = newname + ROp.Getname(ROp) + ")";
  1257.      temp->Nameit( newname );
  1258.      Dispatch->Push( *temp );
  1259.      delete temp;
  1260.      return Dispatch->ReturnMat();
  1261.   }
  1262.  
  1263. VMatrix& Diag( VMatrix& ROp )
  1264.   {  // make a diagonal matrix from a vector or the principal diagonal of
  1265.      // a square matrix, off diagonal elements are zero
  1266.      ROp.Garbage("Diag");
  1267.  
  1268.      if( ROp.r != ROp.c && ROp.c != 1)
  1269.        Dispatch->Nrerror("Diag: ROp is not square or is not a column vector");
  1270.      Dispatch->Inclevel();
  1271.  
  1272.      VMatrix *temp = new VMatrix;
  1273.      *temp = Fill( ROp.r, ROp.r, 0.0);
  1274.  
  1275.      int rr = ROp.r;
  1276.      int cc = ROp.c;
  1277.      for( int i=1; i<=rr; i++)
  1278.         temp->M(i,i) = ( rr==cc ) ? ROp.m(i,i) : ROp.m(i,1);
  1279.  
  1280.      strtype newname = "diag(";
  1281.      newname = newname + ROp.Getname(ROp) + ")";
  1282.      temp->Nameit( newname );
  1283.      Dispatch->Push( *temp );
  1284.      delete temp;
  1285.      return Dispatch->DecReturn();
  1286.   }
  1287.  
  1288. VMatrix& Shape( VMatrix& ROp, int nrows )
  1289.   {   // reshape a matrix into n rows, nrows must divide r*c without a
  1290.       // remainder
  1291.       ROp.Garbage("Shape");
  1292.  
  1293.       long nr = (long) nrows;
  1294.       long lr = (long) ROp.r;
  1295.       long lc = (long) ROp.c;
  1296.  
  1297.       if( lr*lc % nr )
  1298.         Dispatch->Nrerror("Shape: nrows divides r*c with a remainder");
  1299.  
  1300.       Dispatch->Inclevel();
  1301.       VMatrix *temp = new VMatrix;
  1302.       *temp = ROp;
  1303.  
  1304.       long lrc = (lr*lc) / nr;
  1305.       if ( nr >= 32768 || nr < 1 || lrc >= 32768 || lrc < 1 )
  1306.          Dispatch->Nrerror("Shape: reshape produces invalid indices");
  1307.       temp->r = nrows;
  1308.       temp->c = (int) lrc;
  1309.  
  1310.       strtype newname = "shape(";
  1311.       newname = newname + ROp.Getname(ROp) + ")";
  1312.       temp->Nameit( newname );
  1313.       Dispatch->Push( *temp );
  1314.       delete temp;
  1315.       return Dispatch->DecReturn();
  1316.    }
  1317.  
  1318.  
  1319. ///////////////////////////// summation functions // maxs and mins
  1320.  
  1321. //// typedef enum margins { ALL, ROWS, COLUMNS } margins;
  1322.  
  1323. VMatrix& Sum( VMatrix& ROp, margins marg )
  1324.   { // sum(ROp, ALL) returns 1x1
  1325.     // sum(ROp, ROWS) returns 1xc
  1326.     // sum(ROp, COLUMNS) returns cx1
  1327.     ROp.Garbage("Sum");
  1328.  
  1329.     int r,c;
  1330.     double sum;
  1331.     switch( marg ){
  1332.       case ALL: r = 1; c=1; break;
  1333.       case ROWS: r=1; c=ROp.c; break;
  1334.       case COLUMNS: r=ROp.r; c=1; break;
  1335.       default:   Dispatch->Nrerror("Sum: invalid margin specification");
  1336.     }
  1337.     VMatrix *temp = new VMatrix("t",r,c);
  1338.  
  1339.     switch( marg ){
  1340.       case ALL: sum = 0.0;
  1341.                 for( int i=1; i<=ROp.r; i++ )
  1342.                    for (int j=1; j<=ROp.c; j++) sum += ROp.m(i,j);
  1343.                 temp->M(1,1) = sum;
  1344.                 break;
  1345.       case ROWS:for(  j=1; j<=c; j++ ){
  1346.                    sum = 0.0;
  1347.                    for (int i=1; i<=ROp.r; i++) sum += ROp.m(i,j);
  1348.                    temp->M(1,j) = sum;
  1349.                 }
  1350.                 break;
  1351.       case COLUMNS:for(  i=1; i<=r; i++ ){
  1352.                       sum = 0.0;
  1353.                       for (int j=1; j<=ROp.c; j++) sum += ROp.m(i,j);
  1354.                       temp->M(i,1) = sum;
  1355.                     }
  1356.                     break;
  1357.       }
  1358.  
  1359.       strtype newname = "sum(";
  1360.       newname = newname + ROp.Getname(ROp) + ")";
  1361.       temp->Nameit( newname );
  1362.       Dispatch->Push( *temp );
  1363.       delete temp;
  1364.       return Dispatch->ReturnMat();
  1365.   }
  1366.  
  1367. VMatrix& Sumsq( VMatrix& ROp, margins marg )
  1368.   { // sumsq(ROp, ALL) returns 1x1
  1369.     // sumsq(ROp, ROWS) returns 1xc
  1370.     // sumsq(ROp, COLUMNS) returns cx1
  1371.     ROp.Garbage("Sum");
  1372.  
  1373.     int r,c;
  1374.     double sum;
  1375.     switch( marg ){
  1376.       case ALL: r = 1; c=1; break;
  1377.       case ROWS: r=1; c=ROp.c; break;
  1378.       case COLUMNS: r=ROp.r; c=1; break;
  1379.       default:   Dispatch->Nrerror("Sumsq: invalid margin specification");
  1380.     }
  1381.     VMatrix *temp = new VMatrix("t",r,c);
  1382.  
  1383.     switch( marg ){
  1384.       case ALL: sum = 0.0;
  1385.                 for( int i=1; i<=ROp.r; i++ )
  1386.                    for (int j=1; j<=ROp.c; j++) sum += ROp.m(i,j)*ROp.m(i,j);
  1387.                 temp->M(1,1) = sum;
  1388.                 break;
  1389.       case ROWS:for(  j=1; j<=c; j++ ){
  1390.                    sum = 0.0;
  1391.                    for (int i=1; i<=ROp.r; i++) sum += ROp.m(i,j)*ROp.m(i,j);
  1392.                    temp->M(1,j) = sum;
  1393.                 }
  1394.                 break;
  1395.       case COLUMNS:for(  i=1; i<=r; i++ ){
  1396.                       sum = 0.0;
  1397.                       for (int j=1; j<=ROp.c; j++) sum += ROp.m(i,j)*ROp.m(i,j);
  1398.                       temp->M(i,1) = sum;
  1399.                     }
  1400.                     break;
  1401.       }
  1402.  
  1403.       strtype newname = "sumsq(";
  1404.       newname = newname + ROp.Getname(ROp) + ")";
  1405.       temp->Nameit( newname );
  1406.       Dispatch->Push( *temp );
  1407.       delete temp;
  1408.       return Dispatch->ReturnMat();
  1409.   }
  1410.  
  1411. VMatrix& Cusum( VMatrix& ROp)
  1412.   { // cumulative sum along rows
  1413.     ROp.Garbage( "Cusum");
  1414.     VMatrix *temp = new VMatrix("t",ROp.r, ROp.c);
  1415.  
  1416.     double sum = 0.0;
  1417.     for( int i=1; i <= ROp.r; i++){
  1418.        for( int j=1; j <= ROp.c; j++){
  1419.           sum += ROp.m(i,j);
  1420.           temp->M(i,j) = sum;
  1421.        }
  1422.     }
  1423.     strtype newname = "cusum(";
  1424.     newname = newname + ROp.Getname(ROp) + ")";
  1425.     temp->Nameit( newname );
  1426.     Dispatch->Push( *temp );
  1427.     delete temp;
  1428.     return Dispatch->ReturnMat();
  1429.  }
  1430.  
  1431. VMatrix& Mmax( VMatrix& ROp, margins marg )
  1432.   { // Mmax(ROp, ALL) returns 3x1
  1433.     // Mmax(ROp, ROWS) returns 3xc
  1434.     // Mmax(ROp, COLUMNS) returns cx3
  1435.     ROp.Garbage("Sum");
  1436.  
  1437.     int r,c;
  1438.     switch( marg ){
  1439.       case ALL: r = 3; c=1; break;
  1440.       case ROWS: r=3; c=ROp.c; break;
  1441.       case COLUMNS: r=ROp.r; c=3; break;
  1442.       default:   Dispatch->Nrerror("Mmax: invalid margin specification");
  1443.     }
  1444.     VMatrix *temp = new VMatrix("t",r,c);
  1445.  
  1446.     switch( marg ){
  1447.       case ALL: double mmax = ROp.m(1,1);
  1448.                 int maxr = 1, maxc=1;
  1449.                 for( int i=1; i<=ROp.r; i++ )
  1450.                    for (int j=1; j<=ROp.c; j++)
  1451.                        mmax = (mmax < ROp.m(i,j))?
  1452.                 maxr=i, maxc = j, ROp.m(i,j) : mmax;
  1453.                 temp->M(1,1) = (double) maxr;
  1454.                 temp->M(2,1) = (double) maxc;
  1455.                 temp->M(3,1) = mmax;
  1456.                 break;
  1457.       case ROWS:for(  j=1; j<=c; j++ ){
  1458.                    double mmax = ROp.m(1,j);
  1459.                    int maxr = 1; maxc = j;
  1460.                    for ( i=1; i<=ROp.r; i++)
  1461.                        mmax = (mmax < ROp.m(i,j))? maxr = i, ROp.m(i,j):mmax;
  1462.                    temp->M(1,j) = (double) maxr;
  1463.                    temp->M(2,j) = (double) maxc;
  1464.                    temp->M(3,j) = mmax;
  1465.                 }
  1466.                 break;
  1467.       case COLUMNS:for(  i=1; i<=r; i++ ){
  1468.                    double mmax = ROp.m(i,1);
  1469.                    int maxr = i, maxc = 1;
  1470.                    for ( j=1; j<=ROp.c; j++)
  1471.                       mmax = (mmax < ROp.m(i,j))? maxc = j, ROp.m(i,j):mmax;
  1472.                       temp->M(i,1) = (double) maxr;
  1473.                       temp->M(i,2) = (double) maxc;
  1474.                       temp->M(i,3) = mmax;
  1475.                    }
  1476.                    break;
  1477.       }
  1478.  
  1479.       strtype newname = "max(";
  1480.       newname = newname + ROp.Getname(ROp) + ")";
  1481.       temp->Nameit( newname );
  1482.       Dispatch->Push( *temp );
  1483.       delete temp;
  1484.       return Dispatch->ReturnMat();
  1485.   }
  1486.  
  1487.  
  1488.  
  1489. VMatrix& Mmin( VMatrix& ROp, margins marg )
  1490.   { // Mmin(ROp, ALL) returns 3x1
  1491.     // Mmin(ROp, ROWS) returns 3xc
  1492.     // Mmin(ROp, COLUMNS) returns cx3
  1493.     ROp.Garbage("Sum");
  1494.  
  1495.     int r,c;
  1496.     switch( marg ){
  1497.       case ALL: r = 3; c=1; break;
  1498.       case ROWS: r=3; c=ROp.c; break;
  1499.       case COLUMNS: r=ROp.r; c=3; break;
  1500.       default:   Dispatch->Nrerror("Mmin: invalid margin specification");
  1501.     }
  1502.     VMatrix *temp = new VMatrix("t",r,c);
  1503.  
  1504.     switch( marg ){
  1505.       case ALL: double mmin = ROp.m(1,1);
  1506.                 int maxr = 1, maxc = 1;
  1507.                 for( int i=1; i<=ROp.r; i++ )
  1508.                    for (int j=1; j<=ROp.c; j++)
  1509.                        mmin = (mmin > ROp.m(i,j))?
  1510.                            maxr = i, maxc = j, ROp.m(i,j):mmin;
  1511.                 temp->M(1,1) = (double) maxr;
  1512.                 temp->M(2,1) = (double) maxc;
  1513.                 temp->M(3,1) = mmin;
  1514.                 break;
  1515.       case ROWS:for(  j=1; j<=c; j++ ){
  1516.                    double mmin = ROp.m(1,j);
  1517.                    int maxr = 1, maxc = j;
  1518.                    for (int i=1; i<=ROp.r; i++)
  1519.                        mmin = (mmin > ROp.m(i,j))? maxr=i, ROp.m(i,j):mmin;
  1520.                    temp->M(1,j) = (double) maxr;
  1521.                    temp->M(2,j) = (double) maxc;
  1522.                    temp->M(3,j) = mmin;
  1523.                 }
  1524.                 break;
  1525.       case COLUMNS:for(  i=1; i<=r; i++ ){
  1526.                       double mmin = ROp.m(i,1);
  1527.                       int maxr = i, maxc = 1;
  1528.                       for (int j=1; j<=ROp.c; j++)
  1529.                           mmin = (mmin > ROp.m(i,j))? maxc=j, ROp.m(i,j):mmin;
  1530.                       temp->M(i,1) = (double) maxr;
  1531.                       temp->M(i,2) = (double) maxc;
  1532.                       temp->M(i,3) = mmin;
  1533.                     }
  1534.                     break;
  1535.       }
  1536.  
  1537.       strtype newname = "min(";
  1538.       newname = newname + ROp.Getname(ROp) + ")";
  1539.       temp->Nameit( newname );
  1540.       Dispatch->Push( *temp );
  1541.       delete temp;
  1542.       return Dispatch->ReturnMat();
  1543.   }
  1544.  
  1545. //////////////////////////////// elemenatary row and column operations
  1546. ////////////////// note these functions operate directly on ROp
  1547.  
  1548. void Crow( VMatrix& ROp, int row, double c)
  1549.   {  // multiply a row by a constant
  1550.      ROp.Garbage("Crow");
  1551.  
  1552.      if( row < 1 || row > ROp.r )
  1553.         Dispatch->Nrerror("Crow: row index out of range");
  1554.  
  1555.      for( int i=1; i<=ROp.c; i++) ROp.M(row, i) = c*ROp.m(row,i);
  1556.      return;
  1557.    }
  1558.  
  1559. void Srow( VMatrix &ROp, int row1, int row2)
  1560.    {  // swap rows
  1561.       ROp.Garbage("Srow");
  1562.       if( row1 < 1 || row1 > ROp.r || row2 < 1 || row2 > ROp.r)
  1563.         Dispatch->Nrerror("Srow: row index out of range");
  1564.  
  1565.       double tmp=0;
  1566.       for( int i=1; i<=ROp.c; i++){
  1567.           tmp = ROp.m( row1, i);
  1568.           ROp.M(row1, i) = ROp.m(row2, i);
  1569.           ROp.M(row2, i) = tmp;
  1570.       }
  1571.       return;
  1572.    }
  1573.  
  1574. void Lrow( VMatrix &ROp, int row1, int row2, double c)
  1575.   { // add c*r1 to r2
  1576.  
  1577.       ROp.Garbage("Srow");
  1578.       if( row1 < 1 || row1 > ROp.r || row2 < 1 || row2 > ROp.r)
  1579.         Dispatch->Nrerror("Lrow: row index out of range");
  1580.  
  1581.       for( int i=1; i<=ROp.c; i++)
  1582.         ROp.M( row2, i) = ROp.m(row2,i) + c*ROp.m(row1,i);
  1583.  
  1584.       return;
  1585.    }
  1586.  
  1587. void Ccol( VMatrix& ROp, int col, double c)
  1588.   {  // multiply a col by a constant
  1589.      ROp.Garbage("Ccol");
  1590.  
  1591.      if( col < 1 || col > ROp.c )
  1592.         Dispatch->Nrerror("Ccol: col index out of range");
  1593.  
  1594.      for( int i=1; i<=ROp.r; i++) ROp.M(i, col) = c*ROp.m(i,col);
  1595.      return;
  1596.    }
  1597.  
  1598. void Scol( VMatrix &ROp, int col1, int col2)
  1599.    {  // swap cols
  1600.       ROp.Garbage("Scol");
  1601.       if( col1 < 1 || col1 > ROp.c || col2 < 1 || col2 > ROp.c)
  1602.         Dispatch->Nrerror("Scol: col index out of range");
  1603.  
  1604.       double tmp=0;
  1605.       for( int i=1; i<=ROp.r; i++){
  1606.           tmp = ROp.m( i, col1);
  1607.           ROp.M(i, col1) = ROp.m(i, col2);
  1608.           ROp.M(i, col2) = tmp;
  1609.       }
  1610.       return;
  1611.    }
  1612.  
  1613. void Lcol( VMatrix &ROp, int col1, int col2, double c)
  1614.   { // add c*c1 to c2
  1615.  
  1616.       ROp.Garbage("Scol");
  1617.       if( col1 < 1 || col1 > ROp.c || col2 < 1 || col2 > ROp.c)
  1618.         Dispatch->Nrerror("Lcol: col index out of range");
  1619.  
  1620.       for( int i=1; i<=ROp.r; i++)
  1621.         ROp.M( i, col2) = ROp.m(i, col2) + c*ROp.m(i, col1);
  1622.  
  1623.       return;
  1624.    }
  1625.  
  1626.  
  1627.  
  1628.